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Abstract 

In this paper, a new polynomial chaos based framework for analyzing linear systems with probabilistic parameters 
is presented. Stability analysis and synthesis of optimal quadratically stabilizing controllers for such systems are 
presented as convex optimization problems, with exponential mean square stability guarantees. A Monte-Carlo 
approach for analysis and synthesis is also presented, which is used to benchmark the polynomial chaos based 
approach. The computational advantage of the polynomial chaos approach is shown with an example based on the 
design of an optimal EMS-stabilizing controller, for an F-16 aircraft model. 


I. Introduction 

In this paper we study the problem of designing state feedback controllers for linear time invariant systems with 
probabilistic system parameters. Such systems in continuous time can be defined as 

x = A(A.)x + B(A.)u, (1) 

where A S l d is a vector of uncertain parameters, with joint probability density function p( A). Matrices A{ A) € 
r m , b (A) G K nxm are system matrices that depend on A. Consequently, the solution x := x(t, A) G 1™ also 
depends on A. The objective is to design a state-feedback law in the form u = Kx, which stabilizes the system 
in some suitable sense, where K G R m x 71 . Thus, we are looking to obtain a constant deterministic gain K that 
stabilizes the systems in ([]}• with probabilistic uncertainty in A. 

Stability of dynamical systems of the type 

x = A(A.)x, (2) 

have been extensively studied in the framework of stochastic dynamical systems. Depending on the nature of A, 
two approaches are commonly used. If A is Gaussian white noise, the solution process is a diffusion process and is 
analyzed using theory of Markov processes ID, 11 and Ito calculus 0. If A is not Gaussian white noise, and has 
well defined samples properties, the system in 0 can be analyzed by ordinary rules of calculus. Such systems can 
be considered to be a collection of ordinary deterministic differential equations upon which a probability measure 
has been induced by the parameter process A. The second class of systems are often what we actually encounter 
in engineering problems. The Gaussian white noise case is a mathematical abstraction and admits richer set of 
analysis tools. However, this should not be the motivation for assuming white noise process for systems with 
random parameters as it is a non-trivial matter and should be considered carefully 0. 

Early work on stability analysis of linear systems with randomly time-varying parameters was done by Rosen- 
bloom (5j], where he studied first order linear systems with stationary Gaussian parameter process. His work focussed 
on stability properties in terms of the moments. Due to the stationary properties of A, the moments and hence the 
asymptotic properties of system could be studied explicitly. This work was extended by Bertram and Sarachik J6l 
for linear diagonal systems and were first to provide a Lyapunov framework for studying such systems with random 
parameters. The parameters were also considered to be Gaussian and recovered Rosebloom’s result, with a more 
general approach. They also considered the problem where A is random but piecewise constant. At the same time, 
Kats and Krosovskii 0 independently provided a Lyapunov framework for analyzing stability in probability and 
moments, for systems where the A process is a stationary Markov process with finite number of states. Bharucha 
E) extended the work by Kats et al. by exploiting the fact that a stationary Markov parameter process admits a 
piecewise constant property in the linear system. Bharucha also showed that for such systems, asymptotic stability 
implies exponential stability. Palmer E) significantly extended the work by Kats et al. and Bharucha on linear 
systems with Markov coefficients by applying Markov chain theory and exploiting the induced piecewise constant 


property. The piecewise constant property was also exploited in the works of Morozan ESI and Soeda et al. El. 
The above described literature focussed on systems for which the solution can be obtained in closed form. More 
general linear systems of the form x = (Ao + Ai(A))x, where Ao is a stability matrix and Ai(A) is matrix 
whose non-zero coefficients are stationary, ergodic processes with almost surely continuous sample functions, were 
studied by Kozin E2- He provided sufficient conditions for asymptotic stability with probability one, for such 
systems. This work was refined further ESI, E3-E3- For related work when A is Gaussian white noise, please 
refer to ESI, ED and references therein. 

In this paper, we focus on systems where A is a vector of random variables, which is a simpler problem than 
randomly time-varying parameters. The problem of analyzing systems with uncertain, but constant, parameters have 
also been addressed in the robust control literature. In that approach, the support of A is assumed to be polytopic, 
and stability is analyzed for parameter combinations along the vertices of the poly tope fl20i - f22l . This is the so 
called “worst-case” approach for stability analysis. In this paper, we present new stability analysis and control design 
methods, which ensure exponential mean square stability for systems with probabilistic system parameters. These 
are developed in the polynomial chaos framework lf23l . A comparison with worst-case quadratic stability approach 
El, for systems with uniformly distributed A j25l is also performed, which highlights the conservativeness of the 
worst-case approach. 

The paper is organized as follows. First, we present a brief background on polynomial chaos theory and its 
application to linear systems with random parameters. This is followed by the main contribution of this paper, 
captured in propositions 1, 2, and 3; and theorems 1, 2, and 3. For benchmarking, we also present randomized 
algorithms for system analysis and design. We demonstrate the computational efficiency of polynomial chaos 
framework with an example based on an F-16 aircraft model. The paper ends with a summary section. 

II. Polynomial Chaos Theory 

Polynomial chaos is a non-sampling based method to determine evolution of uncertainty in dynamical system 
with probabilistic system parameters ll26l . Very briefly, a general second order process X(u>) £ £ 2 (G, T, P) can 
be expressed by polynomial chaos as 

OO 

X(u) = ^2xi(j)i(A(uj)), (3) 

2—0 

where uj is the random event and 0,(A(uj)) denotes the polynomial chaos basis of degree p in terms of the random 
variables A(w). (fis a probability space, where D, is the sample space, T is the er-algebra of the subsets 
of O, and P is the probability measure. According to Cameron and Martin l27l such an expansion converges in the 
£2 sense for any arbitrary stochastic process with finite second moment. In practice, the infinite series is truncated 
and X(ui) is approximated by 

N 

X(lj) « X(w) = 

2 = 0 

The functions are a family of orthogonal basis in £ 2 ( 12 , J~,P) satisfying the relation 

E [Mi] ■■= [ <MA)^(A)p(A) dA = hfa, (4) 

where d,; 7 is the Rronecker delta, h, is a constant term corresponding to J v (f>jp(A) dA, 72a is the domain of 
the random variable A(w), and p(A) is a probability density function for A. 

Polynomial chaos theory is becoming an useful framework to study control systems with random parameters 

m-m. 


A. Application to Dynamical Systems with Random Parameters 

With respect to the dynamical system defined in ijT}. the solution can be approximated by the polynomial chaos 
expansion as 

N 

x(t, A) « x(t, A) = Y^ Xi(t)<f)i{ A), 

2 = 0 


( 5 ) 


( 6 ) 

(7) 


where the polynomial chaos coefficients Xi £ R n . Define <f> ( A) to be 

$ = $(A) := (<MA) ••• <MA)) T , and 

* n = * n (A) :=*(A)®I n , 

where I n £ R nxra is identity matrix. Also define matrix X £ R nx ( JV + 1 \ with polynomial chaos coefficients Xi, 
as 

X = [x 0 ■■■ x N ] . 

This lets us define x(t, A) as 

x(f,A):=X(f)$(A). (8) 

Noting that x = vec (x), we obtain an alternate form for ©, 

x = vec (x) = vec (X<&) = vec (I n X<&) 


= (3> T <8) J„)vec (X) = &^x pc , 

where x pc := vec (X), and vec (•) is the vectorization operator 031 . 

Since x from © is an approximation, substituting it in (0 we get equation error e, which is given by 

e := x — A(A.)x 

= ^n^pc ~ A(/\)<l> n X pc . 

Best C 2 approximation is obtained by setting 

(e<j>i) := E [efa} =0, for i = 0,1, • • ■ , N. 

Upon simplification, we get a set of n(N + 1) deterministic ordinary differential equations in x pc , 


x pc — E 








Using the following properties of Kronecker product 

A 8 (£ 8 C) = (A ® B) ® C, 
(A 0 B)(C ®D) = (AC) 8 (BD), 

we have the following propositions. 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 


Proposition 1: For <f> n , and $ as defined earlier 

3* 8 ) 8 In- 

Proof: 

$ 8 = 3> 8 (<f> T 8 In) = (<& 8 ^> T ) 8 In 

= ($ ■ 1 ) 8 (1 • $ T ) 

= ($ 81 )( 18 $ T ) 

= ($$ T )8/„ 


5 In 
■ In. 


(14) 


Corollary 1: 


Proof: Straight forward using 0. 


($$ T )8/ n ] = diag (E \4>i]) 8 I n 


(15) 











( 16 ) 


Proposition 2: For any matrix A £ M" x " and $ as defined earlier 


Proof: 


$<g> (A*?) = (W T ) # A. 

* ® (a*?) = (I N+1 *) <g> (a*£) 

= (W <8> 

= (I N+X ** T ) ® (AI n ) 
= ($$ T )<g> A. 


With (fl4l) and (fl6] >. we can write ( fTTl i as 


Xp C — A pc x pc , 


where 


A — F 




and <f>, A depend on A as defined earlier. 

Let us also introduce notations A := E 

E 4> m R4>^ , to compactly present the results below. 


«f> A4> t 


1 E ( 4>4 > T ) 

, B := E 


(17) 


(18) 




, Q := E <&„Q<1>„ , and R := 


III. Stability 

In this section, we study the exponential stability of the second mean for the dynamical system in (0. The 
equilibrium solution is said to possess exponential stability of the m th mean if 3 S > 0 and constants a > 0,3 > 0 
such that ||a:o|| < <5 implies Vi > to 

E [||cc(i; xoM)\\Z\ < /8E [\\x 0 \C] e~ a ^\ (19) 

It can be trivially shown that the dynamical system in 0, with random variables A, is exponentially stable in the 
2 nd mean, or exponentially stable in the mean square sense (EMS-stable), if 3 a Lyapunov function V(x) := x T Px, 
with P = P T > 0, and a > 0 such that 


V 


< —aE [V} . 


( 20 ) 


As the sample properties are well defined, the derivation is obtained by mimicking the proof for deterministic 
systems. Can. (f20l) implies 

clE [V] 


dtt 


< —aE [V] 


Recall for P = P 1 , 


=> E [x T Px] < E [xqPxo] e““ ( ‘“ to) . 
X m in{P)\\x\\l < x T Px < A max (P)||at|| 2 , 


^A min (P)E [||a:||i] < A max (P)E [||* 0 |ll] e- a{t ~ t0 \ 
E [||®||1] < k(P)E [||*o||1] 

EMS-stability condition in the polynomial chaos framework is presented next. 
Proposition 3: For any matrix M £ R m x and 4> rl as defined earlier 

M4> T n =4> T m {I N+1 ®M). 


( 21 ) 

















Proof: 


= * T ® M = (& T I N+ 1) ® (I m M) 

= ($ T ®I m ){I N+1 ®M) 

= <(Iat +1 ®M). 

■ 

Corollary 2: 

$„M t = (Jjv + 1 ® M T )$ m . (22) 

Proof: Take transpose of dTTk ■ 


Theorem 1: The dynamical system in 0 is EMS-stable if 3 P = P T > 0 and a > 0 such that 


A P + VA + q-E 


* (b 1 


P < o, 


where V := I 


N+l 


P. 


Proof: With V(x) := x T Px , and P = P T > 0, 

V = x T Px + x T Px 

= x T ( A t P + PA ) x 


(a t p + pa'j 

= x^ c (& n A T P&„ + ®nPA&l) x pc- 

Using (OH and (l22l >. replace P'i' 1 and <& n P by < P j/P and 'PT*,, respectively. This gives us 


V = x T vc (§> n A T *lV + V^ n A T ^ T N ) x pc 


Similarly, 

V < — aE [V] is equivalent to 


V — x Px — x pc *& n P*& n Xp C — X pc *f* n *f* rl PXp C . 


(23) 


A T V + VA + a E V<0. 


Ean.(l23l> is a convex constraint in P and a (pg. ii, ED). 

In our previous work j30|, we presented stability conditions with Lyapunov function defined as V(x pc ) := 
Xp C Px pc , and studied stability of the deterministic system in (fTTk The dynamics and the Lyapunov function were 
both deterministic and standard Lyapunov arguments were used to analyze stability. Stability of (0 was then inferred 
by showing that 

lim x pc (t) -s- 0 => E [||®|G ->• 0. 

t— koo 

Here we analyze the stability of 0 directly, with Lyapunov functions of the type V(x) := x T Px. 


IV. Controller Synthesis 

In this section, condition for a state feedback controller K that achieves EMS-stability for the system in 0. is 
presented. Let V(x) := x T Px, where P = P T > 0, be the Lyapunov function that certifies this. The closed-loop 
system with control u = Kx is therefore 


x= (A(A) + B(A)K^x. 


(24) 








Theorem 2: Closed loop system in (l24l ) is EMS-stable if 

yA T + Ay + w t b t + bw + a yE 


* 


< 0 , 


(25) 


where W := I N+1 <8> W, y := I N+1 <8 Y, Y := P 1 , and W := KY. 
Proof: With V(x) := x T Px, where P = P T > 0, 


V = x 1 


= x 


pc 


A t P + PA + K t B t P + PBK 




<&r,K T B T P&Z 


<&„PBK$z 


Using m and d22l) . we get 


V = xi 


+ + K 1 SnB 1 $>Vt +V* n B*ilC 


where K. := I 


N+l 


' K. Therefore, E 


V 


< — aE [V] is equivalent to 


A T T + VA + JC t B T V + VBIC + aE 


* 


V < 0. 


(26) 


The above equation is nonconvex in V and JC and can be convexified using the well known substitutions 
Y := P _1 , and W := KY. These substitutions can be written in terms of V,K.,y, and W as 

W = In+i <8 W = I N+1 I N+1 0 KY 
= {In+i <8 K)(In+i 8 Y) 

= icy. 

It is also straightforward to show V = y~ 1 and K, = WJ/ 1 ■ Substituting these in 
by y, we get the following convex inequality in Y.W and a 


and pre-post multiplying 


tT 


yA^ +Ay + W T B T + BW + ayE <0. 


The performance parameter a can be maximized, along with (l25l >. as a generalized eigen-value problem (GEVP) 

(pg. li, ED)- 

Theorem 3: A fixed gain K that minimizes 


( x 1 Qx + u 1 Ru)dt 


Wo 


(27) 


subject to it = Kx and dynamics given by ([]} can be synthesized in the polynomial chaos framework by solving 
the following convex optimization problem 

max try 


subject to 


yA T + Ay + w t b t + bw 
y 
w 


y 

-Q -1 

0 


W T 

0 

R 1 


< 0. 


Proof: The problem is equivalent to solving minp jf trP, such that [x 1 Px\ < E x T (Q + K T RK)x 
In the polynomial chaos framework, the constraint can be written as 

E [$„(A + BK) t P®1 1 + E [*„P(A + BK)$l] + E [*„(Q + K T RK)<&1 1 < 0 




























Using (THT) and ( El ), we get 


A r V + VA + K. t B T V + VBK, + Q + 1C T RK. < 0. 

Substituting Y = P 1 and W = KY , or equivalently y = V" 1 and W = K,y, and pre-post multiplying by 
y we get 

yA T + Ay + \y T B T + bw + yQy + w t rw < o . 

By redefining the cost function in terms of Y and making use of the Schur complement, we arrive at the result. ■ 
In our previous work ll30l . four controller synthesis algorithms were presented, among which was the formulation 
u = Kx, with constant deterministic gain K. The synthesis condition for this case was presented in terms of a 
bilinear matrix inequality, which could not be convexified. Here, we present a convex formulation for the same 
synthesis problem. 

Next we present Monte-Carlo based algorithms for stability analysis and controller synthesis. These will serve as 
benchmarks for the polynomial chaos based algorithms. Let {A s }f =1 ~ p( A) be samples drawn from distribution 
p( A), and the sample trajectories be x(t, A s ). Define X mc := [x(t, Ai) • • • x{t, As)], and x mc := vec (X mc ). 
Recall that for a function F( A), the expected value E [f’( A)] can be approximated as 

E [F(A)] := [ F(A)p(A)dA 

JPa 

1 s 

« ^y'f(A s ), with A s ~p(A). (28) 

^ S=1 

Therefore, EMS-stability condition with Lyapunov function V(x(t, A)) := x T (t, A)Px(t, A), x s := x(t, A s ) 
and A s := A(A s ), can be approximated as 

S= 1 S=1 

^ ( A ^ p+pA s) x s ^ -<*\ Pxs > 

s —1 s =1 

1 s 

- ^ x T s {A t s P + PA S + aP)x s < 0, 

S= 1 

^a;^ c diag (A T S P + PA S + aP) s s=l x mc < 0, 



or 

A t s P + PA s +aP < 0, for s = 1, • • • ,S. (29) 

Similarly, a feasible gain K can be obtained by satisfying 

YA t s + A S Y + W t B t s + B s W + aY < 0, (30) 

for s = 1, ■ ■ • , S, with B s := B(A S ). Finally, the optimal gain for cost function given by (l27l ). can be obtained 
by solving the following optimization problem 

max try 


subject to 


y Ai + A S Y + W 1 B{ +b s w 


w 


Y 

w 


Q 1 0 <0, (31) 

0 -R- 1 

for s = 1, • • • , ,S'. Conditions in ( |30| ) and ( |3TI ) are derived by approximating the expectation integral using ( |28| ) and 
following the argument presented in the derivation of i 




V. Example 

The plant considered here is an F-16 aircraft at high angle of attack [1361 , with states x = [V a q 6 T ] t , 
where V is the velocity, a the angle of attack, q the pitch rate, 6 is the pitch angle, and T is the thrust. The controls, 
u = [5 e Sth] T , are the elevator deflection 5 e , and the throttle 5th- This linear model is obtained about a = 35°, 
where the moderately high angle of attack causes inaccurate modeling of aerodynamic coefficients and thus results 
in parametric uncertainty. The A(A) and B matrices are given by 


' 0.1658 

-13.1013 

-7.2748(1 + 0.4A) 

-32.1739 

0.2780 


0 

-0.0706 ' 

0.0018 

-0.1301 

0.9276(1 + 0.4A) 

0 

-0.0012 


0 

-0.0004 

0 

-0.6436 

-0.4763 

0 

0 

, B = 

0 

-0.0157 

0 

0 

1 

0 

0 


0 

0 

0 

0 

0 

0 

-1 


64.94 

0 


(32) 

Similar to the model in 136j, the terms in parenthesis in the .4(A) matrix are assumed to be uncertain and 
are functions of a single random variable A, uniformly distributed over [—1,1]. This uncertainty is due to the 
uncertainty in the damping term C xq , which is difficult to model in high angle of attack. The uncertainty is 
assumed to be distributed uniformly by ±40% about the nominal values —7.2748 and 0.9276 respectively. Fig.([]} 
shows the variation of ||P*|| obtained by solving the optimization problem outlined above for both Monte-Carlo 
and polynomial chaos based formulation. The objective function is defined by Q := diagf [1 E — 1, 1E2 , IE — 
2, 1E0, IE — 3]), and R := diag([l£’2, 5 E — 3]). In the Monte-Carlo approach, the problem is solved with sample 
sizes 5,10, 50,100, 200, 500,1000,1500, and 2000 respectively. To capture the confidence in the solution for each 
sample size, the optimization problem is solved 100 times with new samples. Fig. © shows the mean ||P*|| and its 
standard deviation. We can observe that the mean ||P* || converges and the confidence in the solution also improves as 
we increase the sample size. Fig. © also shows the convergence of ||P*|| with increasing order of polynomial chaos 
approximation, solved with approximation order 1, • • • ,20. Since the polynomial chaos approach is a deterministic 
framework, there is no variability in the answer for a given approximation order. The x-axis in fig.© is the actual 
number variables in the optimization problem that is solved and is obtained from CVX lf37l /SDPT3 11381 . This is 
proportional to the number of samples in Monte-Carlo approach and the order of approximation in polynomial chaos 
approach. Comparing the two plots we see that the trend with polynomial chaos is consistent with mean trend from 
Monte-Carlo. However, to get a high confidence solution, the Monte-Carlo approach requires to solve a problem 
with an order of magnitude higher in size. With polynomial chaos approach we can get close to the converged 
solution with fairly low order approximation, and thus offers a computational advantage for the class of problems 
considered here. Although, fig.© shows this trend for the specific system, we can expect to see this advantage 
for general systems as well. This is because, in general, polynomial chaos is computationally more efficient than 
Monte-Carlo for characterizing uncertainty in dynamical systems. 

VI. Summary 

In this paper we presented new convex conditions for EMS-stability of linear dynamical systems with probabilis¬ 
tic uncertainty in system parameters. These results were obtained using the polynomial chaos framework and are 
similar to the well known results for deterministic systems. We applied this framework to design EMS-stabilizing 
controllers for an F-16 aircraft model and demonstrated the computational advantage of polynomial chaos based 
approach over Monte-Carlo based approach for analyzing dynamical systems with random parameters. 
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